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Abstract 

We perform a functional expansion of the fidelity between two uni- 
tary matrices in order to find the necessary conditions for the robust 
implementation of a target gate. Comparison of these conditions with 
those obtained from the Magnus expansion and Dyson series shows 
that they are equivalent in first order. By exploiting techniques from 
robust design optimization, we account for issues of experimental fea- 
sibility by introducing an additional criterion to the search for control 
pulses. This search is accomplished by exploring the competition be- 
tween the multiple objectives in the implementation of the NOT gate 
by means of evolutionary multi-objective optimization. 

To appear at J. Phys. A: Math. Theor. 

1 Introduction 



One of the challenges of coherent control of quantum systems is to 
achieve high fidelity in the presence of errors and/or noise that may 
be difficult or impossible to reduce by directly applying more precise 
controls. This situation is also exacerbated for systems with either 
complex underlying interactions or composed of heterogeneous ensem- 
bles. An example is the goal of achieving broadband inversion of spin 
ensembles [U |2] and more generally broadband excitation of spin sys- 
tems. Such needs have lead to the development of composite pulses 
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[31 in O [6] which have been studied in the demanding field of quantum 
computing O O [9] . This technique was extended to include the use 
of shaped pulses [10l[Tl]. A more systematic approach is through the 
application of quantum optimal control [12l [131 [HI [E] > with concep- 
tual foundations lying in the control landscape topology for generating 
unitary transformations |16[ [T7] . The Magnus expansion is commonly 
used t20j to assess the robustness of implementing a unitary gate, in 
contrast to utilizing the Dyson series. The reason for this preference 
seems to arise from the fact that the Magnus expansion maintains uni- 
tarity while the truncated Dyson series does not. However, the fidelity 
between unitary matrices as an objective function is more naturally 
expressed in terms of the Dyson series. The next section validates the 
use of the Dyson series as an appropiate method for the assessment 
of robustness and shows how this relates with the Magnus expansion 
and the series expansion of the fidelity. 



2 Formulating Conditions for Robust- 
ness 



The optimal control of quantum gates for a system of N discrete 
levels may be formulated in terms of the fidelity between two unitary 
operators as a scalar cost function 



where W is the target unitary operator [18j, and the unitary evolution 
operator obeys the Schrodinger equation 



with h absorbed in the Hamiltonian and T being the target time. A 
perturbation in the Hamiltonian H ^ H + SH implies a variation in 
U, which can be assimilated in an auxiliary operator V, defined such 
that 



J(C/(r,0)) ^ -Re{Tr[W^U{T,0)]) 



(1) 




(2) 



u ^u' = uv. 



(3) 



The Schrodinger equation ([2]) implies 



i-V = 5H{t)V, 



(4) 
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with 6H{t) = U^{t, 0)6H U{t, 0). The solution of this equation can be 
expressed in terms of the Dyson series as 

oo 

V{T,0) = re-*/o''5^W'it = 1 + (5) 

n=l 

where Pn are the time-ordered integrals 

rT r-t\ i-tn-i 

Pn = dti dt2... / dtnSH{h)6H{t2)...SH{tr,), (6) 

Jo Jo Jo 

with, for example Pi = dtSH{t). 

Defining AC/(T,0) = U'{T,0) - U{T,0), the following expression 
can be obtained 



AUiT,0) = U{T,0)Y,{-irPn- (7) 

n=l 

Equation ([7]) also can be written as a functional Taylor expansion, 

oo ^ 

AC/(T,0) = ^-5"C/(r,0), (8) 

n=l 

which implies the following identity 

5"[/(r,o) = n!c/(r,o)Hrp„. (9) 

To proceed we define the action of the following brackets that 
specify the Hermitian and anti-Hermitian operators as well as the real 
trace 

{X)h ^ ^(X + Xt) (10) 

{X)a ^ \{X-X^) (11) 

{X)o ^ j^Tr[{X)H] = ^Re{Tr{X)), (12) 
such that we can verify the following identities 

X = {X)a + {X)h (13) 

{{X)a)h = {{X)h)a = (14) 

{{Y)h{X)a)o = (15) 

{Y{X)a)o = {{Y)a{X)a)o (16) 

{Y{X)h)o = {{Y)h{X)h)o (17) 
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With these definitions, the fidehty maybe written as 

j([/(r,o)) = (Tyt[/(T,o))o, (18) 

and the functional Taylor expansion of the fidelity takes the form 

oo ^ 

AJ{U{T,0)) = Y,-{W^^''U{T,0)^ (19) 

n=l 

The first order term becomes 

(w^5UiT,0))^ = (^{W^U{T,0))Ai-i)Pi)^, (20) 

which can be used to define the condition for the regular critical points 
of J{U) [la IT] as 

{W^U{T, 0))a = ^{W^U{T, 0) - U{T, 0)^W) = 0. (21) 

Thus, only the Hermitian part of W^U(T,0) remains at the critical 
points 

W^U{T,0)Uucai = {W^U{T,0))h. (22) 

This implies that the expansion (I19p evaluated at the regular critical 
points becomes 

oo 

Aj(c/(r,o))| 

critical — 

n=2 

which can be used to identify the relevant factors {{—i)^Pn)H that 
depends on the control field. Elimination of the first order term char- 
acterizes a critical point and elimination of higher orders can be used 
as indicators of robustness. The robustness condition extracted from 
the second order term is 

{{-i?P2)^ = 0. (24) 

The Magnus expansion of a unitary operator [23], around the tar- 
get U (T, 0) can be written as 

oo 

W = U{T,0)eMY.'^k), (25) 
k=i 
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where 0^ are Hermitian operators that can be written in terms of P„ 
^i5j according to the identity 

oo oo 

^ ink = log(l + ^i-irPn)., (26) 



fc=l n=l 



such that 



ifli = -iPi (27) 

1^2 = -P2 + \P1 (28) 

= iP3+'^P'l-'^{PiP2 + P2Pi) (29) 

*^^4 = Pi-\{PiP3 + P?,Pi)-\pi+ (30) 

i(PiPiP2 + P1P2P1 + P2P1P1) - \pt (31) 

The criteria for robustness is based on sequential ehmination of Q.^ 
[201 12H [22] , starting from the leading term 

Oi = 0. (32) 

This condition implies Pi = 0, which seems to be unrelated with the 
condition in (f24l) . Applying condition ([32]) the leading terms of ififc 
become 

J7i = (33) 

iOa = (-i)^P2 (34) 

jQa = (-^')P3 (35) 

jQ4 = {-i)^PA + \{-i?Pl (36) 

Extracting the Hermitian part of each term 

^VL2)h = {{-ifP2)H (37) 

(iJ)3)H = {{-ifP-i)H (38) 

(^J^4)h = {{-ifP^)H + \{{-i?Pi)H (39) 

and recalling the anti-Hermiticity of each term of the the Magnus 
series one obtains 

= {{-ifP2)H (40) 

= {{-ifPz)H (41) 

= {{-ifP,)H + \{{-i?Pl)H. (42) 



indicating that J7i = implies {{—i)'^P2)H = 0, showing the complete 
equivalence of conditions (j32p and ()24p . Moreover, r^i = also im- 
plies {{-ifP3)H = and the relation {{-i)'^P4,)H = -\{{-ifPi)H, 
which is useful for writing the two leading terms characterizing the 
robustness according with ()23p as 



Pi = (43) 
{Pi)H = 0. (44) 

However, the last condition can be further simplified considering {P2)h '' 
from (j24p leading to the conditions on the Dyson series 

Pi = (45) 
P2 = 0. (46) 

Moreover, these conditions are consistent if and only if 

ni = (47) 
= 0. (48) 

Convergence of the Dyson series for N-level systems is assured if 
the field is bounded for a finite interaction time T [26]. In contrast, 
convergence of the Magnus expansion demands more severe conditions 
|27j . The convergence is not relevant if the analysis is done in terms of 
the infinitesimal form of the perturbation Hamiltonian 6H. However, 
in practice the perturbation Hamiltonian is finite, implying that the 
Magnus expansion may not necessarily converge. For this reason, any 
proposed robust implementation must be numerically verified for a 
finite range of small perturbations, as performed later in this paper, 
for a specific case. 



3 NOT Gate 

This section is concerned with the implementation of a robust NOT 
gate against off-resonant perturbations based on the lowest order con- 
dition Pi = for robustness. A general form for a single-qubit Hamil- 
tonian with a resonant interaction is 

H = l^oas + ft{t)ai cos{ujot), (49) 
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where ft(t) is the time-dependent modulated Rabi frequency. The 
corresponding time dependent Schrodinger equation is 

Hip = i^V- (50) 

A transformation to the rotating frame wih remove the diagonal term 
of the Hamiltonian. In the present case we consider a more general 
transformation with a time-dependent phase that can be controlled 
using a chirped pulse. The proposed transformation is 

iP = U^ (51) 

with 

U = e~^(^o*-*W)'^3/2, (52) 

where $(t) is the accumulated off-resonant phase generated by chirp- 
ing the pulse. The Schrodinger equation becomes 

{U^HU -iU^—U)^ = i—^. (53) 
The first term is explicitly given as 

U^HU = -iooas + (e^^o*'^^ + g-^^o^'^s )e-^('^o*-*(*))'^3 (54) 

= lujoas + (e**(*)-3 + e(-2i-ot+^*(t))-3 ) (55) 

The last term is highly oscillatory ujq >> $(i). This leads to a 
generalized rotating wave approximation as 

Z^t^Z^^ l^o^3 + M^^e'*W-3. (56) 

Additional comments concerning the conditions specifying this ap- 
proximation are found in|Aj The second term on the left side of (j53p 
becomes 

- ^^^^ = 2 - = (-2^° + 2 ^''^ 

and the final form of the Schrodinger equation in the rotating frame 
is 



2 dt 



U3 + in(t)cTie**W'^=')^I' = z^^, 



(58) 



or, in a more expanded form 



'ld^(t) 1 \ d 

-^^(T3 + 2 ^^(O(cos $ C7i + sin $ aa) 1 * = z^^-, 



(59) 



where both the off resonance phase ^{t) and the Rabi frequency il{t) 
are considered as the control functions. The term = is the 
shift of the resonant frequency as a function of time. It is important to 
note that the adiabatic condition was not required in the formulation 
above. 

The simplest off-resonance perturbation is a constant, which can 
be modeled as 

SH = ^^3, (60) 
taking the following form in the interaction picture 

dH{t) = ^-^u\t)a^U{t). (61) 

Consideration of this constant perturbation is restrictive, but a signif- 
icant degree of robustness remains even for more general off-resonance 
perturbations including the important case of random perturbations 
as shown at the end of this section. 

The NOT gate, up to a global phase, can be generated with the 
following unitary operator that may be implemented with a simple 
square pulse 

NOT = e^''''\e=o^^, (62) 

where the time of interaction occurs on the interval t € [0,T]. Im- 
plementation of the robust gate can be achieved by introducing a 
composite unitary operator, which can be written as 

U{e) = e^^"W(e), (63) 

with V{e) = et'^(^)°'iet^(^)^3, such that the following boundary con- 
ditions are imposed 

L(0) = i?(0) = L(7r) = i?(7r) = 0, (64) 

in order to ensure that V{0) = V{7r) = 1. The associated Hamiltonian 
can be determined as 

H = i%U\ (65) 
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where the energy is measured in units of This Hamiltonian may 
be exphcitly evaluated as 

H = -l{l + L'{e))ai - \ M0 + L{e)]R'{e)a2 - \ cos[0 + L(e)]i?'(0)a3, 

(66) 

which can be compared with (j59p to identify 



n{d) = ^{i + L'{e)Y + sm^[e + L{e)\R'{ef (67) 

v{e) = -cos{e + L{e)R'{9)). (68) 

The boundary conditions (j64p ensure that the target operator wiU 
be achieved, but it is necessary to impose additional conditions in 
order to assure that Q{6) is zero at the boundaries and therefore 
avoid sharp corners at the beginning and at the end of the pulse. The 
additional boundary conditions are 

L'(0) = L'(7r) = -1 (69) 

Furthermore, we also require the modulation of the Rabi frequency 
ft{6) to be symmetric around tt/2. The complete set of boundary 
conditions is satisfied by the following harmonic forms 

n 



k=l 



R{e) = ^bksm{2ke), (71) 
k=l 

for integer n > 1 and Y12=i = ^, where the coefficients are assumed 
to satisfy \ak\ < 2ir and \bk\ < 27r. Dropping the symmetry of Q,{9) 
would result in R{0) = Yl^=i bk sm{k9). 

The minimization of Pi can be carried out through the numerical 
minimization of the following associated objective function for a finite 
number of harmonics n 



•^5J?(°^l'"2, ■■■an,bl,b2, ■■■bn) 



u{ey(j^u{e)de 







(72) 



with an = 1 — X]fc=i ^fc- '^^^ minimization of this cost function ensures 
a robust implementation of the target gate but there is the additional 
desire to implement the pulse with limited experimental resources. 
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This problem can be addressed by introducing additional objective 
functions associated with the shape of ft{9) and I'iO). It is natural 
to consider a Gaussian form as model for Q,{9) due to its smoothness 
and analytical properties 

0(^) = Aoexp(-i^^^). (73) 

The Gaussian function obeys the following differential equation 

The following integral of the residual measures the degree of dissimi- 
larity with respect to a Gaussian, thus, it can be used as a cost function 
for minimization 



(HI, _I^) 



(75) 



The same technique can be applied to the chirp frequency v{0), but in 
this case the simplest chirp available in the laboratory is linear, which 
leads to the minimization of the following cost function 

J, = r de\v"{e)\. (76) 



JO 

In studies on robust design optimization, accounting for practical 
feasibility is typically carried out by introducing an additional crite- 
rion into the search for a control |28l I29j . We choose to follow this 
scheme, and therefore strive to explore the competition between 
and Jq by means of Pareto optimization, employing the MO-CMA-ES 
algorithm (for details, see [HI). A series of runs, considering five con- 
trol parameters (ai, 02, 61, 62, ^3) and a population of 100 candidate 
solutions, produced the Pareto front shown in Figure [TJ The front has 
an interesting shape, revealing a non-trivial conflict between J^f^ and 
Jq, yet allowing a reasonable trade-off, e.g., in the knee point shown 
in Figured]. Setting a threshold of J^fj < 0.0005, the minimum value 
of Jq in the distribution is found to be 3.33, which corresponds to the 
knee point. This point is characterized by the harmonics shown in 
Table [H The corresponding plots of the Rabi modulation and chirp 
frequency are shown in Figures [2] and [3l The robustness of this im- 
plementation is evident in Figure IH which shows the loss of fidelity 



10 



6 




0.00 0.05 0.10 0.15 0.20 0.25 0.30 

Figure 1: The attained approximate Pareto front describing the competition 
between J^^ and Jq. The figure depicts the set of 300 non-dominated points, 
constructed by means of Pareto ranking of 10 fronts obtained in 10 individual 
runs. The enlarged dot is the knee point, chosen as the best compromise. 



n 


an 




1 


0.896833 


3.0578 


2 


0.302287 


0.429276 


3 


-0.207685 


0.0881475 



Table 1: Set of harmonics that produce J^fj = 0.00023 and Jq = 3.33. 

as a function of the perturbation in comparison with an implemen- 
tation with a square pulse (|62p . This figure also shows the response 
of the fidelity to a random Gaussian perturbation (instead of a con- 
stant perturbation) applied along the Pareto front. This figure implies 
that even though most random perturbations result in a fidelity lower 
than that obtained with a constant perturbation, there are some cases 
where the fidelity is actually higher. 

Upon deploying the MO-CMA-ES algorithm on a second bi-criteria 
minimization problem, where the chirp objective competes against 
the coherent average objective J^^j, the redefined goal is observed to 
possess no conflict, which is not a priori evident. This is the case 
because the introduction of an additional objective usually leads to a 
conflict with the original objective(s). It is thus possible to obtain a 
robust implementation by practically ignoring the chirp. A solution 
of this kind is shown in Table [2j The respective plot of the Rabi fre- 
quency modulation is shown in Figure[5]and the response of the fidelity 
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Figure 2: The modulation of the Rabi frequency n{6) for the harmonics 
shown in Table [1] 




Figure 3: The chirp frequency z/(6') for the harmonics shown in Table [H 
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Figure 4: Fidelity plots of the square pulse (dashed lines) and optimized 
implementation of the NOT gate (solid curve), for the harmonics shown in 
Table [TJ The fidelity is shown as a function of the perturbation amplitude 
normalized with respect to the maximum amplitude of Q{6). The dots show 
the response of the robust optimized implementation to a random Gaussian 
perturbation with mean e and standard deviation |e|/2. 
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n 


an 


bn 


1 


2.35701 


5. X 10^^ 


2 


-1.56989 


1. X 10-^ 


3 


0.21289 


1. X 10-s 



Table 2: Set of harmonics that produce Jgfj = 10 ^ and Jq = 0.00002. 




4 2 4 

e 



Figure 5: Modulation of the Rabi frequency fl{6) for the harmonics shown 
in Table El 

to perturbation is practically identical to that in Figure |H Overall, 
this is a good illustration of scenarios where Pareto optimization may 
confirm or dispute the existence of competition between objectives 
where intuition may be misleading. 

In this section we proved that the robustness condition Pi = can 
be used to find a modulation of the Rabi frequency for a qubit system. 
This condition can in principle be applied to higher dimensions but in 
these cases the field cannot be directly extracted from the Schrodinger 
equation and more involved numerical techniques, such as D-MORPH 
[16] . need to be applied. 



4 Conclusions 

We presented an analysis of the functional Taylor expansion of the fi- 
delity between unitary operators with the aim of extracting conditions 
for the robust implementation of target gates. This expansion was 
written in terms of the Dyson expansion and compared with the Mag- 
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nus expansion of unitary operators. The first order term of the fidehty 
expansion is zero when the target is achieved while higher orders are 
associated with the robustness of the implementation. The analysis 
of the Magnus expansion differs because the robustness is associated 
with all the higher order terms including the first one. We showed 
that the second order term of the fidelity expansion and the first or- 
der of the Magnus expansion are equivalent measures of the fidelity 
because eliminating either of them implies elimination of the other. 
This analysis was extended to the next leading order term showing 
additional connection between the fidelity and Magnus expansions. 

Furthermore, we considered the implementation of the NOT gate 
as a case study, while taking into account an additional objective to 
obtain more desirable control pulse shapes. The competition between 
the objectives was successfully identified by means of an evolution- 
ary multi-objective algorithm, allowing for a systematic exploration 
of the objectives and the nature of the conflicts. One case revealed 
an interesting Pareto front, with a promising trade-off area, while the 
competition in the other case was demonstrated to be non-existent 
despite initial expectations of conflict. This case study constitutes an 
example of a scenario where Pareto optimization is needed for balanc- 
ing possibly conflicting gate control objectives, and at the same time 
assessing the validity of initial assumptions, often led by intuition. 



A Generalized Rotating Wave Approx- 
imation 



Section[3]makes use of a form of the rotating wave approximation. The 
justification of this approximation is best understood by analyzing the 
exact case and its limitations. The Schrodinger equation ()53p without 
any approximation is 



cr3-F-fi(t)cos(a;ot)(cos($(t)-a;ot)cri+sin($(t)-cjot)o-2)^' = i^*. 



-{1 + L'{e)) = fl{9)cosiu}o9)cos{<l>{e) -LVoO) (78) 
-sm{e + L{e))R' (6) = n{e)cos{ujoe)sm{^{9) -ujo9) (79) 



1 



2 



(77) 



which can be identified with (1661) to obtain 




(80) 
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such that the modulated Rabi frequency can be extracted as 



= 1 V(l + L'{9)y + sin(^ + mfR'iey. (81) 

I cos(a;oc7j[ 

If the frequency of the free Hamiltonian ojq is on the order of the 
Rabi frequency, there are possible singularities due to the division by 
cos(a;o^)- These singularities can be lifted by proper selection of L{9) 
and R{6), but they are avoided altogether with a shorter pulse (larger 
Rabi frequency). Unfortunately, such a strong pulse may be difficult 
to implement experimentally and in extreme cases the conditions on 
the form of the dipole interaction may not be met. These situations 
are avoided in section 3 by taking a pulse which is weak enough to 
allow evolution under the free Hamiltonian to contain many cycles 
over the control interval. 



B Pareto Optimization with Evolution- 
ary Algorithms 

Pareto optimization aims at simultaneously optimizing a number of 
conflicting objectives, and thereby revealing the Pareto optimal set 
or a region of interest in the trade-off surface between the objectives. 
In this appendix we summarize the principles of Pareto optimization, 
and especially provide some details on our employment of the method. 
Let a vector of objectives in M"^, 

f{x) = {Mx),f2{x),...Jm{x)f, 

be subject to minimization, and let a partial order be defined in the 
following manner. Given any f^^^ G and f^'^^ E M™, we state that 
/(^) strictly Pareto dominates /^^\ which is denoted as 

/ID ^ /I2), 

if and only if the following holds: 

Vi G {1, . . . m} : fl'^ < A 3i G {1, . . . ,m} : ^ < /f ^ (82) 

The crucial claim is that for any compact subset of M™", there exists a 
non-empty set of minimal elements with respect to the partial order 
■< (see, e.g., pO]). Non-dominated points are then defined as the set 
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of minimal elements with respect to the partial order ^. The aim of 
Pareto optimization is thus to obtain the non- dominated set and its 
pre-image in the control space, the so-called Pareto optimal set, also 
referred to as the efficient set. Finally, the Pareto front is defined 
as the set of all points in the objective space that correspond to the 
solutions in the Pareto-optimal set. 

Evolutionary Algorithms (EAs) |3HI32j. are powerful search meth- 
ods, based on natural evolution, which have been successful in treating 
high-dimensional optimization problems. Here, we are especially inter- 
ested in evolutionary multi- objective optimization algorithms (EMOA), 
which have undergone considerable development in the last two decades 
(see, e.g., [33l [Ml ES] ) • Following the broad success of the Covariance 
Matrix Adaptation Evolution Strategy (CMA-ES) (see, e.g., [36]) in 
real- valued single-objective optimization, a multi-objective version has 
been released recently [33 . The Multi-Objective CMA-ES ( MO- 
CMA-ES) is the Pareto optimization approach used in our calcula- 
tions. 

In short, the CMA-ES is an evolution strategy variant that has 
been successful in treating correlations among decision (control) pa- 
rameters by efficiently learning optimal mutation distributions. Ex- 
plicitly, a set of /U search points comprise the evolving population of 
candidate solutions, which correspond to /i independently evolving 
single-parent CMA core strategies. The ultimate goal is thus to ap- 
proximate the Pareto front of the given multi-objective optimization 
problem by means of these fi points. Given the i*'* search point in 
generation (g) of the MO-CMA-ES, x'f\ an offspring is generated by 
means of a Gaussian variation: 

4'^'^ ^M[4'\al'^"c['^) (83) 

The covariance matrices, |C'^^| , are initialized as unit matrices 
and are learned during the course of evolution, based on cumulative in- 
formation of successful past mutations. The step-sizes, li^i^^l > 
updated according to the so-called success rule based step-size control. 
The set of parents and offspring undergoes two evaluation phases, cor- 
responding to two selection criteria: the first criterion is Pareto dom- 
ination ranking, followed by the hypervolume contribution criterion. 
For more details we refer the reader to [37] . 
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